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In this paper we study the hard sphere packing problem in the Hamming space by the cavity 
method. We show that both the replica symmetric and the replica symmetry breaking approxima- 
tions give maximum rates of packing that are asymptotically the same as the lower bound of Gilbert 
and Varshamov. Consistently with known numerical results, the replica symmetric equations also 
suggest a crystalline solution, where for even diameters the spheres are more likely to be found 
in one of the subspaces (even or odd) of the Hamming space. These crystalline packings can be 
generated by a recursive algorithm which finds maximum packings in an ultra-metric space. Finally, 
^ ■ we design a message passing algorithm based on the cavity equations to find dense packings of hard 

' spheres. Known maximum packings are reproduced efficiently in non trivial ranges of dimensions 

, and number of spheres. 
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I. INTRODUCTION 



The problem of packing rigid objects, and spheres in particular, is a fundamental problem which appears across 
disciplines [HQ. In general, the objects could have arbitrary shapes and the ambient space can be an abstract space 
I . A. Given the space and the objects, the main question is that of finding the densest packings. 

In coding theory one is interested in finding an optimal representation of N symbols in binary strings of length 
n, that is A = {0, 1}™ is the Hamming space of dimension n. This optimal coding contains as many as possible 
symbols and ensures that after transmitting through a noisy channel, which at most flips variables, one can 
i ■' recover the original messages. The ratio between the characteristic length of the symbols l c = log 2 N and length of 
the transmitted strings n defines the rate of coding (or packing). The maximum rate of coding is denoted by R. 
Indeed people are interested to know the asymptotic form of R when n, d —> oo and 5 = d/n remains constant. So far 
there is a considerable difference between the best lower and upper bounds for R [3l4l0| . This means that for large n 
1 the best lower and upper bounds for the number of symbols differ by a factor of order 2™. 

Physically, one can consider the set of symbols as a system of identical particles in the Hamming space of dimension 
(f) ' n, interacting by a hard core potential of range d [ll], [Hj]. The aim is then to study the physical states of different 
densities. Clearly for small densities the system is in the liquid phase respecting the translational symmetry. In 
this case it is easy to calculate, for example, the entropy. At higher densities the liquid entropy becomes incorrect 
" (negative) signaling the onset of other stable phases, either crystalline or glassy (T3. 
t-H ■ In this study we formulate the packing problem as a constraint satisfaction problem. We consider N variables (the 
physical particles or the strings of symbols) which take values in A and for each pair of the variables we consider a 
. i-h ■ constraint that forbids overlapping assignments of the two variables. A packing is thus an assignment of the variables 
that satisfies all the constraints. This representation differs substantially with the so called lattice gas models where 
5h , binary variables (representing occupied or empty positions) are defined one each point of the space A. 

The cavity method provides analytical and numerical tools which can be extremely useful in solving optimization 



problems (or constraint satisfaction problems) over random structures [15l - ll9j . In certain cases it is known to provide 
sampling results which cannot be obtained by Monte Carlo Markov chains in subexponential times [e.g. optimization 
problems in the one-step replica symmetry breaking (1RSB) phase]. As an optimization tool it often outperforms 
linear programming methods [20| . The development of such algorithm is, however, by no means obvious due to the 
choice of the representation of the problem and to the need of writing the cavity equations in an algorithmically 
efficient form. 

As we shall discuss in this paper, insights from the application of the cavity method will turn out to be also useful 
in the study of the type of packing problems we are interested in. 

In this paper we study the cavity equations for the packing problem in the replica symmetric (RS) and in the 
one-step replica symmetry breaking (1RSB) approximations [19J. These equations are called belief propagation (BP) 
and survey propagation (SP) equations, respectively [l7l l2lj. In the RS approximation, besides a liquid solution we 
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find a crystalline phase where with higher probability spheres are found in one of the sublattices (even or odd) of the 
Hamming space. This phase has already been observed in Monte Carlo simulations of Ref. [14[. Both the liquid and 
crystalline solutions predict a maximum rate of packing that behaves asymptotically like the best known lower bound 
M. The same result has been obtained in Ref. [l2[ where the liquid entropy is computed in the hypernetted chain 
approximation. To discuss the exactness of the cavity free entropy we also consider some interpolation techniques 
which connect the cavity free entropy to the true entropy of the system [22T - l24l | . 

In the 1RSB case, we provide an approximate solution of the SP equations to calculate the configurational entropy, 
defined as the (log of) the number of pure states or clusters in the solution space. This quantity acquires a nonzero 
value at a clustering transition where the liquid entropy is still positive. The maximum rate of packing that is 
achievable still coincides with the one obtained in the RS approximation. 

Finally we design a message passing algorithm based on the BP equations which allows us to find dense packings in 
not too large dimensions. These packings are typically hard to find by other simple methods like Monte Carlo based 
algorithms. Unfortunately the computation time and memory increase exponentially with the space dimension n, 
making larger dimensions very difficult to explore. To partially overcome this problem we introduce an approximate 
update rule for the message passing algorithm which is restricted to a subspace and makes the computation more 
efficient. This improvement, together with the distributive nature of the algorithm, could help to run the algorithm 
in larger dimensions. We also discuss another iterative algorithm that, given a packing configuration of spheres 
with diameter d, finds another packing for larger diameter d + 1 by increasing the space dimension n. This is a 
polynomial time algorithm generating maximum packings in an ultrametric space. When applied to our problem, we 
find crystalline packings predicted by the RS cavity equations for even sphere diameters. 

The structure of the paper is as follows. In Sec. [IT] we define the problem more precisely and give a summary of 
known results. In Sees. IIIII and llVl we present the BP and SP equations and study their consequences for the hard 
sphere packing problem. In Sec. [V] we study the packing algorithms and Sec. IVII is devoted to extension to the 
g-ary Hamming spaces. Finally the concluding remarks are given in Sec. I VIII In the first two appendixes we give the 
details of calculations for checking the stability of the BP solutions and deriving the SP equations. The interpolation 
methods and some of their properties are presented in Appendix [Cj 



II. DEFINITIONS AND KNOWN RESULTS 



Consider N hard spheres of diameter d indexed by i — 1, ... ,7V and a Hamming space of dimension n. The set 
of points in this space are denoted by A = {0, 1}™ with size |A| = 2™. We index the points in this space by a. A 
point can be represented by a binary vector of n elements € {0, 1}. The Hamming space can be partitioned into two 
subspaces: even and odd. Points in the even subspace have an even number of l's (even parity) and those in the odd 
subspace have an odd number of l's (odd parity). The Hamming distance D(<r,a') between two points a and a' is 
equal to the number of different elements in the binary representation of the two points. For each point a we define 
the set Vd(a) as 
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V d (a) = W\D(a,a') < d}, V d = \V d (a)\ =£(?)• (1) 

i=o ^ ' 

The aim is to find a non-overlapping configuration of spheres such that D(<7i,<jj) > d, for any two spheres i and 
j. The above problem is a constraint satisfaction problem with N(N — l)/2 constraints to satisfy. We index these 
constraints by (ij). A configuration a = {<Ji\i = 1, . . . , N} that satisfies all the constraints is called a solution of the 
problem. The partition function Z counts the number of such solutions, 

/. Y^\\1,,\o,.g^. (2) 

q_ i<j 

where ij, (er^ , a 3 ■) is an indicator function for constraint (ij); it is 1 if D(ai,<r.j) > d and otherwise. The maximum 
rate of packing is defined as 

R= lim -log 2 (iV mo a.), (3) 

n,d— too 71 

with S = d/n = const and N max is the maximum number of spheres such that Z > 0. 

Let us here mention some known lower and upper bounds for R. The Gilbert and Varshamov (GV) lower bound 
0, [j| states that 

2 n 

N>max 2^ TT~ ; (4) 
Vd 
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FIG. 1: Comparing the known lower and upper bounds for the maximum rate of packing R. 



resulting in the following lower bound for the maximum rate of packing 

R > R GV = 1 - H(8), 

where 

ff(*) = -*]og a *-(i-a)io g2 (i-a). 



(5) 



(6) 



Notice that 1 — H (5) is also the Shannon rate for a binary symmetric channel with error probability S. A better lower 
bound is obtained with graph theoretical methods [l(| and gives 



R > R JV = 1 - H(S) + 



log 2 [cnH{5)] 



(7) 



where c is a constant. As far as we know, this is the best lower bound reported for the maximum rate of packing. 
However, it still behaves asymptotically like R GV . The authors in Ref. [12| use the liquid entropy of the system of 
hard spheres to find a maximum rate of packing that is very close to the above lower bound 
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PZ 



1 - H(S) 



log 2 [(21n2)nif(J)] 



On the other side, we have the Hamming upper bound Q 



Vd 



resulting in 



R < R H = 1 - H I - 



(8) 



(9) 



(10) 



The best upper bound for R is obtained by the linear programming methods 0. An overestimate of this linear 
programming bound is 

R < R MRRW = H{^- y/6(l - 5fj . (11) 

In Fig. [T]we have compared the above bounds to show the large gap between the best lower and upper bounds. 
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III. REPLICA SYMMETRIC SOLUTIONS: BP EQUATIONS 



The basic ob jects in the cavity method are the cavity marginals or messages passed along the edges of the interaction 
graph [I?], [Hj]. I R our problem the messages are 2™-component vectors with elements € {0, 1}. There are two 
kinds of messages: (i) The cavity bias Uj_>.j represents the warning that variable i sends to j; uf^j = 1(0) says that 
point cr G A is (not) forbidden by i for variable j. (ii) The cavity field h,-_^ is sum of warnings that i receives in 
absence of j; hf_^,j gives the number unsatisfied constraints, in the absence of j, if variable i takes state a. The cavity 
biases are 

u„, e { e> e A), 4 = { k in) 

In the RS framework we assume that all the packings or solutions belong to the same cluster of solutions in the 
configuration space. Suppose the interaction graph is a tree and cr* is a solution of the problem in the absence of 
variable j. Then we represent the cavity bias by e ff «. Notice that according to our definitions a cavity solution 
uniquely determines the cavity message Uj_>j. The histogram of cavity biases among the solutions is given by 

o, ^ ,,;./„..,.• (13) 

cr 

Assuming a tree interaction graph one can write equations governing the cavity probabilities T/f^j : 



II (£ J ** fa ). 

keV(i)\j \ cr k / 



^i" II [Tjik^i^kH^], (14) 

kev(i)\j 

where V(i) denotes the set of variables interacting with variable i. The above equations are called BP equations 
In the Bethe approximation the entropy density of the system is written as 



1 

S= N 



E As *-E Ai 



(15) 



where Asi and As^ are the entropy shifts by adding variable i and interaction (ij), respectively. For these quantities 
we have 



eAs * = E II [E^fa^AiU^ 

a jeV(i) \ cr' I 



(16) 



e AslJ = E 7 «fa (7 ')^^ = Zij. (17) 



cr.cr' 



In words, Zi is the probability that variable i can occupy at least one point in the Hamming space and is the 
probability that interaction (ij) is satisfied. 

Using the above equations for hard spheres we obtain 



Hkev(i)\j ( 1 Ecr'eVrfO^fe-H, . 
VUj = _ _ ^7 = T^T> ( 18 ) 



Ylcr' TlkeV(i)\j (l J2cr"£V d (a') Vk- 

and 

^ = E II M- E Vfj], ZiJ = l~ E (19) 

c jeV(i) \ cr'eV d (a) I cr,cr' :D{cr,cr')<d 



Now we can try different solutions to the BP equations. 
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FIG. 2: Comparing N%££ (BP) with some lower bounds (LB) from [27] . Up to n — 15 all the lower bounds are exact. 



1. Liquid solution 

The liquid solution is obtained by taking J/jL^ = T) = for any i and j. Evaluating Zi and Zij at the liquid 
solution we can write the BP entropy 

s = l n ( 2 ») + ^_Iln(l-^), (20) 



where Vd = This entropy vanishes at 

n bpl = ! _ 

ln(l - v d ) 



N BPL = 1 - 21n ^ 2 ") . (21) 

max i„/i \ V^ 1 - 1 -/ 



For large n it gives 

A^ x S d ~(21n2)n. (22) 

In Fig. [5] we compare this quantity with some known exact results and lower bounds in small dimensions. The 
maximum rate of packing predicted by the liquid solution of the BP equation is 

^^l-g W + lQg2[(21n2H . (23) 

n 

where we have used 

v d ~ 2- n[1 ~ H{s)] . (24) 
We check the stability of the liquid solution in Appendix [3] we find no continuous glass transition as long as 

21n2<— , (25) 

which is the case for 6 < | as Vd is exponentially small in this region. 

2. Crystalline solution 



In some situations the spheres may prefer one of the sublattices of the Hamming space to the other [12J. Indeed 
for even values of d, a sphere at the origin forbids more points from the odd sublattice than the even one. For high 
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densities the neighboring spheres will be at distance d from the origin, i.e., they occupy points that again belong to the 
even sublattice. The above observation suggests that we could have ordered states in which nearly all the spheres are 
in the even or odd sector of the Hamming space. Here we will consider this situation by assigning different messages 
■q e and rf to points in the even and odd sublattices, respectively. Again we use the fact that all the points in one 
sublattice are equivalent. Now the BP equations are 

(^ — V e n e — V°>n°) N ~ 2 

rf = \L lAH V d^ > (26) 

' 2 n ~ 1 [(l - V<fr] e -Vpf) N - 2 + (1-V£r)° -V°jf) N ~ 2 y y ' 

= (i - vjrf - v^°) N - 2 

' 2™- 1 [(l - V e d if - V°t]°) n - 2 + (1 - Vjif - V°T] e ) N - 2 ) ' 
where V d and V d are given by 

2=0 v / z=o v / 

Notice that rf > and rf = cannot be a solution of the above equations unless for d = n. Indeed to have such a 
solution we need r/ e = -^k=r = yz, which holds only when d = n. We will look for other solutions where vf > 0, i]° > 0. 
In this case we have 

rf (1 - V° A rf - Vi^ f- 2 = r?°(l - Vjrf - V° d rf f-\ (28) 

Let us introduce the rescaled variables v d = V d /2 n ~ 1 , v d = V d /2 n ~ 1 , p e = 2 n ~ l rf and p° = 2 n ~ 1 rf . Using the 
normalization condition p e + p° = 1, we obtain an equation for p e , 

f [l ~ v° dP e - - p e )] N - 2 = (1 - p e )[l - v e dP e - „S(1 - p e )] N - 2 . (29) 

For large n and d a nontrivial solution appears at Nq PC ; see Fig. |3] At this point, the two sides of Eq. l29lhave the 
same slope at p e = 1/2. Thus we find 

< PC = 2 + J^. (30 ) 

The entropy of the crystalline phase is obtained given Zi and Z^ , 

= 2"- 1 [(1 - v\p e - v° dP °) N ~ l + (1 - vfr° - vlp-f- 1 ] , (31) 
Zii = 1 - P e (v e d p e + v° d p°) - p°(v e d p° + v° dP e ). (32) 
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FIG. 4: Comparing the liquid and crystal entropies predicted by the BP equations. 



To get a simple expression for the entropy we do the following approximation: p e = 1 and p° = 0. Notice that when d 
is even, v° d > v e d and from Eq. [29] we see that for large N, p° should be much smaller than p e . In this approximation 
the entropy reads 



This entropy vanishes at 



For large n the leading term is 



fl « ln(2"- 1 ) + ^- -ln(l-«3). (33) 



_21n(2"-i) 



ln(l - v%) ' 



N^v^(2ln2)n. (35) 



For n,d — > oo the ratio Vd/v% approaches to a constant and we recover asymptotically the bound provided with the 
liquid solution. In Fig. 3] we have compared the liquid and crystal entropies for a given n and even d. Later in Sec. 
IV Bl we will introduce an iterative algorithm constructing the above crystalline packings. 



IV. ONE-STEP RSB SOLUTION: SP EQUATIONS 

In the 1RSB framework we assume that there are an exponentially large number Af c ~ of clusters of solutions. 
A cluster of solutions consists of packing solutions in the configuration space that are connected to each other by 
paths of finite Hamming distances in the thermodynamic limit. In each cluster we could have frozen and unfrozen 
variables. Suppose the interaction graph is a tree and consider all cavity solutions (in absence of variable j) that 
belong to a given cluster; in a tree graph, fixing the boundary variables is equivalent to fixing the cluster of solutions. 
It may happen that in all the solutions of a cluster, variable i takes only one state, say a. Then the survey bias Ui^j 
will be e CT and we say variable i is a completely frozen variable in that cluster. We could have partially frozen variables 
that are frozen on a subset V? of the Hamming space. In this case we represent the survey bias by — ^2 cvf e ° 

ignoring the degeneracy of each state. In general we could write Ui^j = J2 a ev J w i->j e °- where > is the 

number of times that variable i appears in state a. A variable that takes all possible values is called an unfrozen 
variable and its survey bias is represented by 0. 

The survey biases defined above depend on the cluster of solutions and change from one cluster to another. For an 
edge (ij) of the interaction graph, we define the histogram of the survey biases among the clusters 

Q^(W)=^>,o + 2 E E (36) 

m— 1 (7i<---<(7 m 
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Notice that compared with the RS case here we have the extra option 0. The survey fields are determined by the 
survey biases: 



H 



Uj= E mii ? u-^V)}- (37) 



k€V( t )\j a ' eV - 



In words, Hf-tj is the minimum number of unsatisfied constraints, in the absence of j, if sphere i takes position a. 
Here we need only to know if Hf^j is zero or greater than zero. So we change the field's definition to 



Ji, E min {i-tikM}} 



HUj = mm I 1, £ min {1-I ik (<r, *')}}■ (38) 
[ kev(z)\j cr 't 

The aim is to write an equation for Jj^j-"' ™. The value of J$_^ is determined by the normalization condition. Variable 
i sends the survey bias Ui^.j = e CT1 + ■ • • + e am to variable j if its survey field is 

UUn = { ?' tf ^ e ^—'^>i (39) 
I 1, otherwise, v 7 

that is if the only choices for variable i are the states in V* = {<7i, . . . , <r m }. We write this probability as 

1 - Prob |J UUj = 1 -OR. (J UUj = ) , (40) 

where [J a "Hf_^ = is the event that at least one of the Hf^j is zero. We have also the condition that variable i 
dose not receive contradictory biases. It means that the cavity field Hi-^-j should not be equal to 1 = (1, ... , 1). This 
probability is given by Prob (\J a V.f_^j — 0) ; so we obtain 



i - Prob ( u aevf nUi = i -or U aeA \ Vf %Ui = o 



^ Prob (U CT ^ =0) ■ 1 > 

One can use the standard tools in the probability theory to find the above probabilities in terms of 77's. The result is 
the so called SP equation and is derived with more details in Appendix [B] Here we simplify the analysis by assuming 
that all partially frozen surveys are zero, that is, only completely frozen and unfrozen surveys are taken into account. 
This would be a reasonable approximation when N approaches N max . We also expect to obtain a larger complexity 
within this approximation. In this ansatz the SP equation reads 



(-i) m E 

li—tj 7 "1 

like 



(42) 



where Vu(o"i, . . . , cr m ) is the union volume of Vd(cri), Vd{<T2), ■ ■ ■ and Vd(er m ). 

In our problem all the edges of the interaction graph are equivalent. Moreover, due to the translational symmetry, 
the surveys r)f_+j do not depend on a. Considering these simplifications, the SP equation can be rewritten in a more 
compact form as 

n= i E^g(y)(i-e a m 

where 

2" 2" 

G (^) = E (-l) m+1 9™(V), G'(V) = E (-ir +1 m<? m (V). (44) 

rn—l m—1 

Here g m (V) is the total number of configurations that m distinct points can take in the Hamming space such that 
the union volume V\j is equal to V. More precisely we have 

9m(V) = E 5 V,Vu{<ri,...,a m )- (45) 



9 



To obtain more explicit results we consider a naive approximation where V u (ai, . . . ,a m ) ~ mVd- This is a good 
approximation for small d/n. We also approximate (1 — raV^f]) N ^ 2 by e _ ( Ar ~ 2 *' mV < 1 ' , to compensate the error made by 
overestimating the union volume Vu(ai, . . . , a m ) for large to. Then the SP equation reads 



On / on 

^^_,(-l) m 1 I ] , -2)\;,,, 



e 

. (46) 



Summing over m we get 



7/ = 1 _ (1 _ e -(JV-2)V d r,)2" ■ ( 47 ) 

Let us see where the above equation admits a nontrivial solution r\ > 0. To this end we will consider the following 

scalings: (N — 2)Vd/2 n = cn and 2 n rj = 1 — g(n, c). For n-Joowe expect to have c — > const and g — > 0. Replacing 
these in the above equation and expanding for small 2"e~ cn< - 1-9 ' we obtain 

g~2 n e~ cn{1 ~ 9) . (48) 

To have a solution for g we need c > In 2 + d (In n)/n with 

(nln2 + c'lnn)e- c ' ln " = 1. (49) 

Increasing N, the frozen variables appear for the first time at Njj? , where 

c '~ 1 + (-L), c ~ln2+ — +o(-), (50) 
\ In n J n \n J 

and therefore 

N^ p v d ~ nln2 + Inn + const. (51) 
This gives the clustering transition, when the solution space splits into an exponentially large number of clusters. 



3. Complexity 

We can compute the complexity in the Bethe approximation as 



1 

N 



£ AS,- AS, 



(52) 



As before, AE^ and AEy are the complexity shifts by adding variable i and interaction (ij). The complexity is 
zero for small densities of the particles where we expect to have a single cluster of solutions. As we increase N, we 
encounter the clustering transition at N c , where E becomes nonzero. The complexity is a decreasing function of N and 
finally vanishes at N max . Here we shall focus on the most numerous clusters, which are not necessarily the relevant 
ones. More accurate results can be obtained by a large deviation study of the problem, which involves computing the 
complexity of different clusters (26j . 

To compute the complexity we need the two quantities AE^ and AEy . Again for the sake of simplicity we only 
work with the completely frozen surveys. In this approximation we get 



AS; _ 



E(-D m+1 E II i- E - 

m=l ci<---<a m jeV(i) \ er'(EV u (eri,...,<7 m ) 



and the link contribution is 



AS, 



= !- E = z k 

a,cr' •.D(a,<j')<d 



(53) 



(54) 
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Using our notation in the previous subsection we can write a more compact form of the complexity for uniform surveys 

E = In 



^Givxi-vn)"- 1 

L v 



?L-±\n[l-2 n V d if]. (55) 



In the naive approximation, where we approximate Vy by mV d , we find 

E = ln[l - (1 - e-^- 1 ^") 2 "] - ln[l - 2 n V d r, 2 ]. (56) 

The complexity vanishes when 

1 - (1 - e -(N-i)v d vf l = (1 _ 2 n V d7? 2 )^ i . (57) 
Using again the scalings (N — l)V d /2 n = cn and 2 n rj = 1 — g(c, n), the above equation can be rewritten as 

I _ e -x = g-icnCl-s)^ a ._ e npa2-c(l-fl)]_ ( 58 ) 

Expanding the left side for the exponentially small a; we obtain 

21n2 = c(l-.g 2 ). (59) 

Notice that g approaches to zero for large n as 1/n, so the above equation suggests 

N^ p ax v d = (2\n2)n + 0(l/n). (60) 

Using Eqs. [47] and [56] we can find the complexity in the naive approximation. In Fig. [5] we have compared this 
complexity with the BP entropy for n = 25 and d = 2, where the naive approximation is expected to work. We see the 
jump in the complexity that happens at the clustering transition ■ Moreover, the complexity is always smaller 
than the RS entropy, as it should be; it is the sum E + s c i uster that gives the total entropy of the system. Here s c i us t er 
is the internal entropy of the clusters. 

Let us summarize the main approximations we used in the above calculations: First, the SP equation ignores the 
soft part of the messages and works only with the hard fields. We know that this could give an underestimate of the 
clustering transition but works well in predicting the satisfiability-unsatisfiability (SAT-UNSAT) transition in some 
constraint satisfaction problems. Next, we neglected the partially frozen states and considered only the completely 
frozen and unfrozen parts to write a simpler expression for the SP equation. This approximation seems reasonable 
as for maximum packings the spheres should be strongly localized by their neighborhood. And finally, we resorted to 
another approximation by replacing the union volume of m distinct spheres with mVd- One can improve on this by 
using an annealed approximation of Mj, 

V u (a 1 ,...,a m )^2 n [l-(l-v d ) m ], (61) 
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where (1 — Vd) m is the probability that a point in the Hamming space is outside the m spheres centered at a%, . . ., 
a m . Further, we can approximate g m (V) by 

g m (V) - ( ^ ) ( v ) I 1 - (! - ^TK 1 " ^T]"- (62) 

The annealed volume is approximately given by mVd when m <C fr-. Moreover, the completely frozen approximation 
we used, suggests that the main contribution in the SP equation comes from small values of m, see Appendix [B] 
Therefore, we do not expect to observe exponentially large deviations in N^ p and N~^ x beyond the naive approxi- 
mation. 



V. SOME ALGORITHMS FOR THE PACKING PROBLEM 



A. Hard spheres in an ultrametric space 

Exact solutions are always useful in that we obtain some insights about more complex problems [28|, [29|. As a 
simple example which can be treated exactly, we consider the packing problem in an ultrametric space. 

An ultrametric space is a metric space where for any three points ii,t2, and 13 we have D{i\ 1 i2) < 
max{Z?(i!, i 3 ), Z?(i 2 , 13)}. Consider a rooted binary tree T n with n generations. The number of points at genera- 
tion I is 2 . The space of points is given by the 2 n leaves of this tree. The distance D(ii,t2) between two points i\ 
and %2 is defined as the number of generations that one should go up in the tree to find the first common ancestor. 

The problem is to find a packing of N hard spheres such that for any two spheres D(ii,i 2 ) > d. The packing 
problem is trivial for d — 1, so we can start from this simple case and try to find packings for larger d. Let us start 
with a binary tree of n\ generations and 2™ 1 leaves, T ni . We put one sphere at each point of this ultrametric space to 
obtain a packing of N = 2 ni spheres with diameter d — 1; see Fig. |6] This is the densest configuration of hard spheres 
for the above parameters. The idea is to increase d by 1 and add the minimum number of necessary generations to 
find a new packing with the new diameter for the spheres. Indeed in an ultrametric space we need only one additional 
generation to remove all the overlaps in the previous configuration. After adding the new generation we are free to 
put each sphere in one of its two descendants. One can continue the above process to find a valid packing for an 
arbitrary value of d. It turns out that these are the densest configurations of the spheres with the parameters n and 
d] the packing indeed partitions the whole space into regions occupied by the spheres. Therefore, given n and d, the 
number of packings and maximum number of spheres read 

/ an — c£ \ 

Z = N\( Z N J(2 d ) JV , N max =T- d . (63) 

The maximum rate of packing will be 

rUM = \og 2 N max = i s 
n 

It is interesting that this is exactly the Shannon rate for a binary erasure channel with error probability 5. 
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FIG. 7: Comparing number of spheres in the crystalline packings generated by the iterative algorithm with some lower bounds 
(LB) from [27j. Up to n = 15 all the lower bounds are exact. 



The above example is one of the exactly solvable packing problems that one can compare the exact N max with the 
BP one computed at the liquid solution: 



« L = 1 - frff^-n) - (2 ^ 2)n2-". (65) 



Here N%££ > N max but we find asymptotically R BPL = 1-5. 



B. Hard spheres in the Hamming space 

Here we can use the same strategy as above to find packings of hard spheres in the Hamming space. The points of 
an n-dimensional Hamming space can be represented by the leaves of a binary tree T„. A configuration of N spheres, 
a = {(ii € A\i = 1, . . . , N}, is a packing of hard spheres with diameter d if it satisfies the following set of constraints 

CMeI^^^^/j), (66) 

where D(<Ji,<jj) is the Hamming distance of points Oi and <7j. We define the energy E[a\ = X)i<j[l — ^i( cr ii <7 j)] as 
the number of unsatisfied constraints. The conflict graph G(a) represents a graph of N nodes where edge (ij) is 
present if the corresponding constraint is not satisfied. To find a packing we do the following steps: 

• We start with N = 2 ni spheres occupying all the leaves of T ni . This is the densest configuration of spheres that 
satisfies C(l). We use q} to represent this configuration of spheres. 

• For t = 2, . . . , d, we start with (<r* -1 , T nt _ 1 ) and find (<r*, T„ t ) such that all the constraints in C(t) are satisfied. 
Set n = nt— i and An t = 0, then 

— If E > 0, add a new generation to T„, that is n = n + 1 and An t = An t + 1. 

— For each sphere find the best value of (Ji(n) 6 {0, 1} and call the new configuration (j 4-1 ^" 4 . If E = 
return ct* = o- t ~ 1,A ™ t and n t = n. 

Notice that every time we add a new generation we have to solve an optimization problem to find the best con- 
figuration, s = {<7i(n)\i = 1, . . . , N}. Given t and An t , we need to find the ground state of the following energy 
function 

m= E <w ( 67 ) 
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FIG. 8: Fraction of locally frozen spheres in the packings generated by the iterative algorithm. 

The aim is to assign different values to the new components Sj and Sj when edge (ij) belongs to the conflict graph. 
In our numerical simulations we use simulated annealing to solve the above optimization problem. 

We observe that for even t the change in the number of generations An t is 1; i.e., the conflict graph is a bipartite 
graph connecting only spheres in different sublattices, odd and even. Consequently, for even d all spheres are in 
the same sublattice of the Hamming space, as expected from the crystalline phase of Sec. MI 21 In Fig. [7] we have 
compared the packings constructed by the above algorithm with the known lower bounds. We also checked the fraction 
of locally frozen spheres in these packings. A locally frozen sphere is one that, fixing the other spheres, cannot be 
locally displaced without violating some of the constraints. Figure [8] displays this quantity as a function of d. We 
see that for even d a considerable fraction of the spheres are frozen. This indicates that finding such packings is very 
difficult by a random sequential adding algorithm [301 ] where spheres are added randomly one by one. 

C. An algorithm based on the BP equations 

Let us again write the general BP equations for the packing problem: 



Wi-tj 



kev(i)\j \ <j k ) 



(68) 



To solve these equations one starts with random initial values for the cavity messages and updates the messages 
iteratively to reach a fixed point. At the fixed point we have the local marginals given by 



jev(») 



(69) 



One can find a packing configuration by decimation, fixing the spheres one by one according to the above local 
marginals; after each decimation one has to run again the BP equations to obtained the new marginals. Here we take 
another approach, by using the local marginals to converge the equations to a polarized fixed point defining a packing 
configuration. We found it useful here to add also some additional random potential Wi{ai) to break the high level of 
symmetry present in the problem. The modified BP equations read 



keV(i)\j \ <? k ) 



(70) 
(71) 
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where cavity messages are biased toward their local marginals. These are the reinforced BP (rBP) equations [3l| and 
the positive parameter r is called the reinforcement parameter. The next step is to write the above rBP equations in 
the limit j3 —} oo, assuming the messages scale like r]i^j{ai) oc e £ m *-»j(°'i). The resulting reinforced max-sum equations 
are 

mi_,.j(o-i) = Wi((Ti) +rm,i{ai) + V max m fc ^(cr fc ), (72) 

uaTTi-w .°-fc : - f ifc(°"^ cr fc)=i 
kev(z)\] 

mf(tri) = Wi{(Ji) + rm,i(ui) + V] max rrij^i(aj). (73) 

We solve these equations by iteration starting from random initial messages and zero reinforcement. The reinforce- 
ment parameter r is increased gradually as r(t + 1) = r(t) + Sr. At the end, the states that maximize the rrii{ai) 
define a packing configuration. 

The time complexity of this algorithm grows as [N(2 n - V d )] 2 . We checked the algorithm by finding some packings 
given in [27]]: for instance, maximum packings with parameters (n = 10, d = 4, N = 40), (n = 11, d = 3, N = 144), 
(ri = 11, d = 5, N = 24), (n = 15, d = 7, N = 32). The main difficulty with large dimensions and number of spheres 
is the computation time and memory. One way to reduce both time and memory is to work with restricted search 
spaces. For example, we may take an initially random search space Aj of size S = 2 n " for each sphere. Given these 
search spaces we run the reinforced max-sum equations to obtain the cavity and local messages. Then for each sphere 
we replace a part of its search space with new states. The above steps are repeated to find better search spaces and 
maybe a packing. 

More precisely, the algorithm starts by initially randomly selected search spaces Aj = {<J^\ ■ ■ ■ , cr> S '}. Then we do 
the following steps: 

• Run the reinforced max-sum equations: The messages are updated for sufficiently large number of iterations T 
using the restricted search spaces Aj. For example, we use T = 100 if we increase the reinforcement parameter 
by Sr = 10~ 2 . 

• Update the search spaces: Compute the local messages {m,i{a)\a € A^} and sort them in a decreasing order to 
put the good states at the beginning of the list. Replace a number of states at the end of this list with other 
states close to the best one. In our simulations we replace half of the search spaces. 

In this way we could, for example, find the maximum packings for (n = 10, d = 4, N — 40) and (n = 11, d = 3, N — 
144) with a search space of dimension hq = 3 and no = 5, respectively. This also allowed us to find some packings 
in larger dimensions, for instance (n = 53, d = 29, N = 10), (n = 54, d = 29, N = 12). These are the best known 
packings (maybe maximum) with those parameters. And in some cases we find denser packings with larger sphere 
diameters: (n = 48, d = 32, N = 4), (n = 51, d = 30, N = 6). 



VI. HARD SPHERE PACKING IN THE g-ARY HAMMING SPACE 



The above results can easily be extended to the y-ary Hamming spaces A = {0, . . . , q — 1}™. There are q n points in 
this space represented by vectors of n elements in {0, . . . , q — 1}. As before, the Hamming distance D(a, a') gives the 
number of different elements in two vectors a and a' . Therefore the number of points at distance less than d from a 
point is given by 



d-i 

n 



V d = \V d (a)\=J2[ n l \ (q-l) 1 . (74) 

In the asymptotic limit n, d — > oo and q,S — d/n finite, we get v d = ^ ~ e~ n ^ 1 ~ Hl >( s ^ where 

H g (6)=Slog q (q-l)-6log q (6)-(l-6)log q (l-6), < 5 < 1 - -. (75) 

In this section we shall only discuss the results obtained within the replica symmetric approximation. The BP 
equations remain the same as in Eq. [HI so we write directly the BP entropy computed at the liquid solution, i.e., 

s=^-^]n(l-v d )+]n(q n ). (76) 
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This entropy vanishes at 



BPL 21n(g") 

max Hi- Vd y [ ' 



which gives a rate of packing that is not asymptotically different from the Gilbert- Varshamov (GV) one, 

R ~l-H q (S) + — . (78) 

n 

But, we know that for q-axy alphabets and q large enough, the rate of packing could asymptotically exceed the GV 



bound |32j. These packings, known as algebraic-geometry codes in coding theory, result in the following lower bound 



when q is a square: 



R > R TV = 1 - S — , 0<6<l — . (79) 



And if q > 49 there exists always an interval in which R TV > R GV . We believe that these are some crystalline or 
ordered solutions that only happen for large values of q. The fact that at S = the rate is smaller than 1 indicates 
that the above packings are restricted to a subspace of A which is exponentially smaller than q n . Moreover, comparing 
R TV with the rate of packing in an ultra-metric space suggest that this subspace is effectively ultrametric; see Sec. 
IV At More precisely, an ultra-metric subspace of size ~ q™ 1 " 1 V?- 1 ' where each sphere occupies ~ q d points would 
result to the same rate of packing as R TV . We remind that in ultrametric spaces both the BP and GV bounds are 
asymptotically exact. Notice that also for q — !• oo one obtains R BPL — R GV = 1 — 5. 

Indeed the number of points that their Hamming distance from a given point is equal to their ultrametric distance 
is U n = ^E§[(<Z — l) n — 1], which is exponentially large as long as q > 3. This is an upper bound for the size of an 
ultrametric subspace. Moreover, given diameter d, we need only a subspace that is ultra-metric up to distances less 
than d. Let us define fi as the probability that a point belongs to an ultrametric subspace. A naive estimation of fi 
can be obtained by 

H~(l-n) Vd - Ud . (80) 

Note that Ud is exponentially smaller than Vd, thus asymptotically fi ~ q- nH t{°) . In this case, each sphere will occupy 
only O(l) points of the subspace, so we recover the GV rate of packing. The above argument says the gap between 
the typical and maximal ultra-metric subspaces is very large, assuming the TV bound comes from packings in an 
optimal ultrametric subspace. 



VII. CONCLUSION 



In summary, we have written the BP and SP equations to study the hard sphere packing problem in the Hamming 
space. Within these approximations we obtained a maximum rate of packing that is asymptotically the same as the 
lower bound of Gilbert and Varshamov. In the RS approximation we also found a crystalline phase where for even 
values of d the spheres prefer to be in one sector of the Hamming space. The BP solutions were stable with respect to 
continuous glass transitions as long as 6 < 1/2. This suggests that phase transitions, if any, would be of discontinues 
type. 

We have also introduced two new algorithms. First a message passing algorithm based on the BP equations which 
finds dense packings of hard spheres in finite dimensions. An approximate scheme is used to reduce the time and 
memory complexity of the algorithm, which can still be improved and hopefully lead to new packing results. An 
almost identical algorithm can be used in continuous spaces. As a proof of concept we used it to find some known 
local dense packings in two-dimensional Euclidean space. 

Second, we introduced an iterative algorithm to find packings of hard spheres starting from small diameters and 
dimensions. For even diameters the algorithm generates packings with all spheres in one subspace of the Hamming 
space, as expected from the crystalline solution of the BP equations. 

There are still some points that need more effort to be clarified. It is of extreme relevance to establish the relation 
between the exact maximum rate of packing and the one provided by the BP equations. This means, for instance, 
that we need to study some interpolating functions between the exact and the BP entropy. Preliminary results are 
discussed in Appendix [Cj We expect the BP bound to be asymptotically exact in the binary Hamming space. 

Finally, it would be nice if one could obtain the algebraic-geometry lower bound for the maximum rate of packing 
in the g-ary Hamming spaces with some physical arguments. These are probably crystalline solutions that cannot be 
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captured within the RSB formalism as the RS entropy is expected to be larger than that of the glassy solutions in 
the RSB phase. 

As a concluding remark, we should mention that the 1RSB study presented here is not complete; the SP equations 
only consider the hard or frozen part of the cavity messages. A more complete study of the 1RSB approximation also 
takes the soft part of the messages into account and asks for the stability of these solutions [33l435| . 
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Appendix A: Stability of the liquid solution 

The liquid solution of the BP equations would be stable as long as the ferromagnetic (linear) and spin glass 
(nonlinear) susceptibilities are finite [33U35| . These conditions can be expressed in terms of the maximum eigenvalue 
of the response matrix M: 



M <y,<r' — 5— -1 3 — { 2 v d ,vfw-,, r ;c a (Al) 



driUj _ f ~^r, ifa'eV d (a); 
drj£n I 2 H $- Vd) , otherwise. 

that has been evaluated at the liquid solution. Then the stability conditions read 

Af PL |A ma:c | - 1, Ni PL \\ max \ 2 = 1, (A2) 

where X m ax is the maximum eigenvalue (in absolute value) of M. 

More precisely, for N < Nf PL and N < Nf PL we would have no continuous phase transition to an ordered and 
spin glass phase, respectively. Notice that these conditions do not exclude discontinuous phase transitions. 

The response matrix M is symmetric with real eigenvalues and orthogonal eigenvectors. Using the translational 
symmetry by even vectors we write the eigenvectors as 

u% = e" k - CT , (A3) 

where a is the binary vector representing a point in the Hamming space and k is a binary wave vector. The eigenvalues 
are obtained by plugging the above expression in the eigenvalue equation 

A < = ~ ^ E "A + 2 nJ d _ V A E ( A4 ) 
o'eV d (a) V d > *>eA\V d (a) 

Then for the eigenvalues we obtain 

A=-— V J***"' -°) + - Vd V e iTk.(<r'-<r)_ / A5) 

2" ^ 2" (2" - Vd) ^ 

<y'eV d (a) X a> a>£A\V d (a) 

The above expression is independent of a and for simplicity we choose a = (0, 0, . . . , 0). Consider the maximum wave 
vector k = (1, !,...,!). The corresponding eigenvalue is 



Obviously, the maximum eigenvalue is bounded by 2vd- Let us approximate the maximum eigenvalue by the largest 
contribution at / = n/2, 

|A m az| — — , 77-T ( n ) ■ (A7) 



2»(2»-V d ) V f 
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Notice that for large n this is still asymptotically equivalent to the trivial bound 2v d . Using the Sterling approximation 
for large n we get 

lXmaxl - 7^(2»-V d )' (A8) 



which according to Eq. IA2I gives 

N* PL v d * - v d ), Nf PL v d * JL(1 - v d f. (A9) 

2 Av d 

Appendix B: Survey propagation equation 

Let us start from Eq. [JT] and find a more suitable form of the SP equation, 

ai a _ 1 - Prob ([j^vf KU, = 1 -OR. U CT£A \y/ KUj = o) 

~ " Prob(a^=0) ' (B1) 

where Vf = {<7i, . . . , er m }. Using the inclusion-exclusion theorem we write 

Prob ( (J = ) = E Prob = °) - E Prob ( fl ^ = ) 

\ (T / (Ti o-i<CT2 \a , =o , i,o , a / 

+ E Prob ( fl ^ =0) -••• + (-l) 2 "- 1 Prob(f|^ J =0) . (B2) 

eri<CT2<er 3 \cr=(Ti,CT2.cr 3 / \ cr / 



Here Prob (fl^,...^ = °J is the probability of having ftf^. = 0, H^j = 0, ... , W^. = 0. 

The probability in the numerator can be rewritten in the same way as 

Prob I |J UUi = 1 -OR- IJ = = I Prob ( U K+i = 1 J + E Prob (k^ = °) I 

- I E Prob f U Wl+i = 1 - AND ' = I + E Prob ( fl "■ v (l ) I • • • (B3) 

In the above equation we have probabilities like 

Prob |J HUj = 1 -AND. p| HUj = 01, (B4) 
\<rey/ creA\v/ / 

which, using the normalization condition, can be rewritten as 

Prob I p| HUj = I - Prob I f] %U 3 = AND. p| HUj =01. (B5) 

\o-eA\W / \a£l^ cr£A\y^ / 

Plugging these into our expression for the numerator, Eq. IB31 we find 

Prob |J UUi = 1 - 0R - U U U = I = I Prob ( U n Uj = 1 ] + E Prob = 0) 1 

I e pr ° b = °) - e pr ° b ( n ^ 



7 -0 AND. 7^ = 

+ e pr ° b ( n ^=o)| 



(B6) 
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Notice that the last term in the first bracket is canceled with the first term in the second bracket. This cancellation 
indeed happens for any two subsequent brackets. Simplifying the above expression we find 

1 - Prob I |J HUj = 1 -OR- (J HUj = 

\aeVf a£A\Vf 

= Prob n k«=o)- e pr ° b ( n = ° and - = ° 

+ J2 Prob ( fl H Uj = -AND. f) HUj=0) -••• + (-l) 2n - |W| Prob(f|^. = 0j . (B7) 

a 1 <a- 2 eA\Vf \aeVf tr=o\ ,<r 2 J \ <J ) 

Now we should write a more explicit relation for Prob (f| a 7if_^j = 01. First note that 



Prob (H°^ = 0) = J] ( 1 " E E ^ € VnW, . . . , o> m ,)\rii^^' I , (B8) 

keV(i)\j \ m'=l cr' x <---<a' m , J 

where Vn( cr i; • • • % a ' m ') is the intersection of sets Vd(ci), ■ ■ ■ , Vd(a' m ,). This is the probability that no variable in V(i)\j 
sends a survey that forbids state o\ for variable i. Here 1(C) is an indicator function for condition C. And in general 

/ \ / 2" n{ffi,...,tr m } \ 

p»m n «^=o = n h-E e . ( bq ) 

\a=a 1 ,...,a m J keV(i)\j \ m' = l o^<-<er' , / 



where 



n{CTi,...,<r m } 

E = E nK-.^ln^,...,^,)^^. (bio) 

o-^<-"<cr' , <tJ<---<<t' , 



Using the above probabilities in Eqs. IB2l and lB7| we obtain the SP equation, 

Em'=0 '(-1)™ Z)<j;<---<<t' ,£A\V/ rifceV(i)\j ~ 7fc->i(°l) • • • ) <T m'>°'l) • • • >°m)) 



°U r ° m = ___ _ -1' nni^lT — — v,w _ _ _ (Bn) 



with 

2" n{ffi,...,tr' ,} 



7^^,...,^,)^ e E fc" CTm "- ( Bl2 ) 

■m" — 1 <t"<" •<(?■"„ 



Appendix C: Interpolating between the BP and exact entropies 

In any approximation it is important to know how well the method approximates the correct behavior. Considering 
the BP approximation, we know that it is exact at least on tree interaction graphs. On loopy graphs, the BP 
approximation may overestimate or underestimate the correct entropy. This in general depends on the nature of the 
interactions and configuration space. 



1. A static interpolation 



The liquid solution of the BP equations is indeed equivalent to treating all the constraints independently and 
replacing a check /^(o-j,^) with its average value when both <Tj and aj are free. To model this solution, on each 
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interaction edge (ij) we introduce auxiliary variables (sij,Sji) besides the original variables. The new variables will 
serve as independent copies of variables (<Tj,<7j) on edge (ij). Clearly if we want to recover the exact partition 
function we should force all the copies to be the same as their originals. Summing all together we write the following 
interpolating partition function: 

m - Ell ( arhw £ hM^ 3i )^^ )+D{ ^ s A , (cd 

a_ i<j \ V Si^Sji J 

with an inverse temperature /? scaling as 1/N. Then we have Z BPL = Z(0) and Z = Z(oo). The above interpolating 
function can be defined for any constraint satisfaction problem that admits a uniform liquid solution. 

Notice that without the normalization factor 1/(1 + e ~P) 2n wc would obtain a replicated partition function Z r ((3) 
that is a convex and decreasing function of /?. In fact Z r (f3) > Z r (oo) = Z, providing a simple upper pound for 
the true partition function for any /?. To get a good upper bound one needs to solve the replicated problem for the 
largest possible (5. Adding the normalization factor destroys the upper bound property but it is necessary if we want 
to recover Z BPL at f3 = 0. 

Let us define variables Xi as 

X i = Y,D(a i ,s ij ). (C2) 



Taking the mth derivative of ln(Z) with respect to (3 we obtain 

d m lnZ _ d m \nZ d m lnZ 

d(-f3) m ~ d(-(i) m ~ d{-/3) m ' 

with 



(C3) 



Z = ^ e -^ x - UhM^sji), Z = Y j e~^ x K (C4) 



One can easily check that at = and for m < 6 

d m In Z _d m In Z a 
d(-P) m ~ d(-P) m ' 



(C5) 



Indeed, it is at in = 6 that for the first time we encounter closed loops connecting three spheres in a diagrammatic 
expansion of the partition function. So for in < 6 we have 

r) m In 7 

frat'« = o- < C6 > 

which means the interpolating partition function is nearly flat at /3 = 0. This is also true for j3 — oo where all distances 
D(ai,Sij) should be zero. 
For m — 6 we have 

d 6 lnZ ^ dHnZ 3 

where Z 3 is the interpolating partition function of three spheres: 



E }Z ^ 2 (rs-Ji)e~ Mr2 

r 3 >dr 1 ,r 2 



(C8) 



Here Si is the number of points at distance I from a given point and Qi ± ,i 2 (l) is the number of points at distance l\ 
from oi and distance h from (T2 when D(a\,a2) — I- More precisely, we have 

Qi,h(l)=( )( V (C9) 
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FIG. 9: The interpolating partition function for three spheres. 

And finally R ri ,r 2 ( r 3i is the number of pairs (s\, S2) that satisfy the following conditions: D(s\,a\) = r x , D(s 2 , 02) — 
r2, and D(s\, S2) = r 3 given D{a\, 02) = I- This number can be written as 

(r). (C10) 

r 

In Fig. [9] we have plotted the typical behavior of Z3 (/3) for some value of n and d. Despite the fact that Z%(0) is 
decreasing with /? we find that the first nonzero term in the Taylor expansion is positive, signaling the nonperturbative 
nature of Z 3 (/3) with respect to fi. 




2. A dynamic interpolation 

Here we use an idea previously introduced in Refs. (23l . l2~i| to replace the interactions in the original problem with 
the effective ones coming from the BP approximation. 

Let us label the edges of an interaction graph by t = 1, . . . , M . Define Eq as the empty set and £ t — {e\, . . . , e t } as 
the set of the first t edges. Now we introduce the following sequence of partition functions 

z t=En*' ff i«)' ( cn ) 

g_ ee£t 

Clearly, the original partition function is given by Z = Zm and 

Z BPL — Zq n*=i M^ti a it ~ a jt))o where the average is taken with respect to the uniform and independent distribution 
of <Xj t and Cj t . Suppose we add edge et+\ connecting nodes (i t +i, jt+i)- Then 

Z t+1 = Z t <J t+1 ) (l + A 4+1 ), A t+1 = (/t+1 , )t ~f t+l) ° , (C12) 

where (I t +i)t = Pt{It+i) is the probability of satisfying constraint It+i when the interaction set is given by £f And 
(7t+i)o = Po(h+i) is the same probability when the interaction set is £q- Moreover, A t+ i = when the interaction 
graph £ t +i is a tree. We rewrite the final entropy as 

M 

logZ = logZ BPL + ^log(l + A t ). (C13) 
t=i 

For hard spheres A t may have different signs depending on the interaction graph. For example, suppose the 
interaction graph is a chain of size L and we add the interaction between the end points. The probability of finding 
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FIG. 10: Absolute value of Al in a chain of repulsive interactions. Filled and empty symbols show positive and negative values, 
respectively. 



the end points at distance r can be obtained in a recursive way as 



1 



2" - V d 



Qry\r') 



(C14) 



'>d 



using our expression for Q h ,i 2 (l) in Eq. [C9j In this case, having Al is enough to know if the BP entropy is an over- 
or underestimation. Figure [TU] displays this quantity for different values of L. We observe that A^ is always negative 
for small d, but for large d and small L its sign alternates as expected for antiferromagnetic interactions. 

Looking again at the BP entropy shows that in general the correction in Eq. IC13I is relevant only if divided by N it 
is still exponential in dimension n. The corrections are irrelevant for instance if A t 's behave as independent random 
numbers with a symmetric distribution and A t = 0(1). This is what we expect to happen in high dimensions. 

It is instructive here to compare the above hard sphere problem with a system of attractive interactions where the 
constraints Jjj (cj, <jj) are satisfied if D{ai, <jj) < d. Figure [TT1 displays A^ in a chain of attractive interactions; A^ is 
positive for small L but approaches an exponentially small negative value for large L and d. In this case, we expect 
to have asymptotically A t > 0. This is similar to the Griffiths inequalities [36] for ferromagnetic systems. Actually, 
using the Fortuin-Kasteleyn-Ginibre (FKG) inequalities [37| one can prove the above statement for soft attractive 

interactions, where /y(<7j,crj) = e~~ Dlyai > a ^ and (3 > 0: 



Consider the following measure on configurations g_ g A where A = {0, 1} 

e£S t 



(C15) 



This measure satisfies the conditions of the FKG inequality; A N is a partially ordered set and ^(s.) is convex, that is 
given two configurations a and a' 



(C16) 



where g_ max = max(<7, a') and g_ mm = min((7, a'). These are bitwise max and min functions. One can use induction 
to show that D{a^ ax ,af ax ) + D{a r ^ in ,af m ) < D(a i ,a j ) + £>(>*, <r£) which proves the convexity of fi(a) for the soft 
attractive interactions. 

Having these conditions the FKG inequality says that for decreasing (increasing) functions f(g_) and g(g_) 



(fg) - (f)(g) > Q, 



(C17) 



where the averages are taken with respect to the measure /i(cr). A function f(a) is decreasing if /(a) > f(cr') for any 
a' > <7. The latter is a bit-wise inequality. 
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FIG. 11: Absolute value of Al in a chain of attractive interactions. Filled and empty symbols show positive and negative 
values, respectively. 



Now consider decreasing functions / = I[D(<Ti, 0) < ZJ and g = I[D(<jj, 0) < lj] where is the zero member of A. 
According to the FKG inequality we have 



(I[D(cn,0) < 0}I[D(a 3 ,0) < I}) > (I[D(ai,0) < 0]) (I[D(aj , 0) < /]), 



(C18) 



which means P t [D(<Ti,<jj) < 1} > Po[D(<Ti, ov) < I] for any i,j, and I. In other words, it is more likely to find two 
spheres closer to each other compared to the uniform measure. As a result, for the soft attractive interactions BP 
provides a lower bound for the log-partition function. 

Actually for the soft interactions = e ± "N- D ( cr i>°j) the partition function reads 



N 



M=0 



N 
M 



,±/3M(l-M/Nhn 



(C19) 



One can exactly solve the problem in the thermodynamic limit and compare the results with the BP predictions to 
confirm the above statement for the attractive interactions. In the case of repulsive interactions the BP and exact 
results asymptotically coincide. 
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